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1.  Introduction 


The  recent  success  of  characteristic-based  upwind  differencing  on  structured  meshes  has  spawned 
significant  research  into  adopting  such  methods  for  use  in  unstructured  adaptive  mesh  flow 
solvers  [1-4].  Several  basic  approaches  for  the  construction  of  characteristic  based  unstructured 
flow  solvers  have  evolved  in  recent  years.  To  first  order,  the  support  stench  employed  in  most  of 
these  methods  are  similar  in  that  they  rely  only  on  next-neighbor  information  in  a  manner  which 
mimics  first  order  structured  schemes.  The  essential  difference  in  the  design  of  these  methods 
lies  in  their  extension  to  higher  order  spatial  accuracy. 

One  approach  towards  obtaining  higher  order  begins  with  constructing  what  is  essentially 
the  structured  difference  stencil  on  the  unstructured  mesh.  References  [2,4],  for  example,  use 
either  Harten’s  modified  flux  approach  [5]  or  van  Leer’s  MUSCL  [6]  reconstruction  to  obtain 
higher  order  accuracy.  Such  schemes  typically  make  use  of  one  additional  point  beyond  the 
nearest  neighbor  in  their  difference  stencils.  Passing  this  information  through  the  unstructured 
mesh  is  the  essential  challenge  faced  in  their  design  and  implementation.  This  noncompactness 
introduces  additional  complexity  into  the  methods  and  usually  requires  additional  storage  to 
overcome.  A  main  benefit  of  this  approach  is  that  the  resulting  schemes  generally  retain  the 
favorable  convergence,  robustness,  and  accuracy  properties  associated  with  upwind  schemes  on 
structured  meshes  [2]. 

A  second  technique  was  introduced  by  Ref.  [1]  using  a  centered  approximation  for  estimation 
of  the  slope  within  each  control  volume,  and  thus  obtains  higher  order  accuracy  without  informa¬ 
tional  inquiries  beyond  the  nearest  neighbor.  The  Green’s  theorem  approach  adopted  by  Barth’s 
linear  reconstruction  [1]  essentially  degenerates  to  the  Fromm  scheme  in  one-dimension  [7]. 
Variations  on  this  theme  have  been  proposed  by  Ref.  [3]  and  others.  Additionally,  the  concept  of 
using  a  nearest-neighbor  based  slope  estimation  was  immediately  generalized  to  use  least  squares, 
constrained  least  squares,  and  higher-order  procedures  [8-10]. 
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Relying  only  upon  next-neighbor  information  maintains  a  compactness  comparable  to  central 
difference  based  schemes  and  permits  the  use  of  edge-based  —  and  other  compact  —  data 
structures.  A  main  strength  of  this  class  of  methods,  is  that  they  avoid  the  ambiguities  associated 
with  the  mapping  of  structured  mesh  based  difference  stencils  onto  an  unstructured  mesh  and  data 
structure  [11]. 

A  second  important  strength  stems  from  the  divestment  of  the  gradient  calculation  from 
the  stencil  associated  with  the  Riemann  solver.  The  method  requires  only  a  best  estimate  of 
the  solution  gradient  within  each  control  volume,  and  while  the  surface  integral  used  in  linear 
reconstruction  still  suffers  on  poor  quality  meshes,  more  elaborate  gradient  estimations  may  not 
Such  procedures  may  provide  dramatically  more  reliable  gradient  estimates  and  suggests  the 
possibility  of  constructing  schemes  with  much  greater  tolerance  for  poor  quality  meshes.  This 
point  becomes  increasingly  important  when  considering  viscous  simulations  [12].  In  such  cases, 
accurate  evaluations  of  the  second  derivatives  often  place  severe  restriction  on  the  quality  of 
acceptable  meshes,  and  the  high  aspect  ratio  cells  necessary  to  resolve  the  stiff  physics  of  high 
Reynolds  number  viscous  flow  may  actually  hinder  accurate  gradient  estimates.  One  goal  of  the 
present  work  is  to  examine  the  scheme  behavior  in  detail  on  the  right  triangular  meshes  often 
employed  to  resolve  viscous  layers. 

An  abundance  of  research  in  the  recent  literature  has  attempted  to  capitalize  on  the  promises  of 
various  reconstruction  techniques.  Nevertheless,  outstanding  questions  concerning  convergence, 
accuracy,  robustness,  and  monotonicity  remain.  The  present  effort  focuses  on  the  simulation  of 
inviscid  flow  with  Roe’s  approximate  Riemann  solver.  Time  integration  is  achieved  via  a  three- 
stage  modified  Runge-Kutta  method.  The  investigations  are  specifically  designed  to  examine 
the  accuracy,  stability  and  expense  of  the  various  reconstruction  algorithms  and  comparisons 
with  structured  MUSCL  schemes  are  provided  wherever  possible.  The  methods  are  examined  on 
meshes  consisting  of  quadrilateral,  right  triangular,  and  equilateral  triangular  elements  to  examine 
the  effects  of  polygon  shape  on  convergence  and  accuracy.  Both  Green-Gauss  and  least  squares 
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reconstructions  are  considered  using  Barth’s  original  Limiter  [1],  Venkatakrishnan’s  limiter  [7] 
and  a  new  directional  limiter  which  is  presented  in  Section  2  of  this  paper. 

The  investigations  are  designed  to  independently  examine  the  methods  in  a  controlled  setting. 
In  addition  to  several  well  known  numerical  test  cases,  these  examples  also  include  problems 
for  which  there  exist  closed  form  solutions  to  the  2-D  compressible  Euler  equations.  Such  cases 
facilitate  quantitative  statements  about  scheme  performance  and  order  of  accuracy. 
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2.  Theoretical  Model 


The  governing  equations  are  written  in  flux-integral  form: 

JJADdA-ifMs  (2» 

where  U  is  the  state  vector  of  conserved  variables  and  F  is  the  tensor  of  flux  density  containing 
the  inviscid  components  of  the  Navier-Stokes  equations.  This  equation  may  be  discretized  on  an 
unstructured  mesh  containing  polygons  of  arbitrary  degree.  Figure  1  shows  a  small  region  of  a 
typical  mesh  with  triangular  elements.  The  present  work  adopts  a  node-based  approach  and  thus 
the  flow  variables  reside  at  the  mesh  vertices.  Hie  solution  procedure  for  Equation  2.1  consists  of 
reconstruction,  flux  quadrature  and  updating  of  the  state  variables.  The  flux  quadrature  requires 
evaluation  of  a  divergence  operator  while  the  reconstruction  step  necessitates  estimation  of  a 
gradient. 


2.1  Formulation 


We  begin  by  deriving  the  edge  formulas  for  the  reconstruction  as  in  Ref.  [13].  This  is  useful 
since  the  analysis  contained  in  subsequent  sections  will  examine  many  of  these  expressions  and 
assumptions  in  significant  detail.  In  Figure  1,  the  perimeter  of  the  polygon  formed  by  the  set  of 
triangles  containing  Vo  as  a  common  vertex  forms  the  closed  path  of  integration  5.  The  polygon 
enclosed  by  this  boundary  has  an  area.  A,  equal  to  the  sum  of  the  areas  of  the  figures  which  share 
the  common  vertex  Vo.  Thus,  the  integral  path  for  any  given  point  is  chosen  to  be  that  described 
by  joining  all  adjacent  points  as  defined  by  the  edge  to  edge  connections.  Applying  Green’s 
theorem  in  the  plane  then  provides  an  evaluation  of  the  gradient  of  a  passive  scalar  <f>  at  the  vertex 


Vo, 


I  Ia  -  j>s  <f>ndS 


(2.2) 
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where  n  denotes  the  outward  facing  local  unit  normal  to  5.  On  the  multifaceted  control  volume 
shown  in  Figure  1,  the  trapezoidal  rule  provides  a  discrete  analog  to  Equation  2.2  and  is  exact  if 
<j>  is  linear.  The  contribution  of  each  edge  K  Vj  to  the  gradient  of  <j>  at  Vo  is  then: 


\  (<Pv,  +  <t>v3)  nv,v, 


where  i,j  vary  cyclically,  i,je{  1 , 2, . . . ,  6}  and  nv.v,  denotes  the  surface  normal  vector  to  edge 
Following  Ref.  [13],  separation  of  the  cell  based  formula  in  Equation  2.3  into  contributions 
from  each  of  the  edges  incident  upon  Vo  permits  reexpressing  the  integration  on  an  edge-basis,  and 
so,  the  trapezoidal  evaluation  of  the  path  integral  in  Equation  2.2  may  be  equivalently  expressed 
in  terms  of  contributions  associated  with  each  node.  For  a  vertex  V}  on  the  perimeter  of  the  cell 
surrounding  vertex  14,  this  term  is 


(nv.Vj  +  nVjvk) 


i,k  are  vertices  adjacent  to  j 


(2.4) 


To  further  facilitate  application  in  an  edge-based  setting,  the  vectors  in  Equation  2.4  may  be 
written  in  terms  of  edges  of  the  centroid  dual  without  changing  the  actual  path  of  integration  S  in 
Equation  2.2.  Ref.  [8]  points  out 


+  nv3vk  =  3ncoyc(b-,k  (2.5) 

where  C,j*  denotes  the  centroid  of  the  triangle  formed  by  vertices  Vi,Vj,Vk.  The  quadrature  of 
Equation  2.2  is  evaluated  with  Equations  2.3  through  2.5  to  obtain: 

f  \/<j>dA  =  <f  <t>nd.S  =  YL  ^v-xncooCo,*  (2-6) 

JA  Js  jev,,...,v6  z 

with  i,  k  once  again  denoting  the  vertices  adjacent  to  j.  Hence,  the  gradient  at  Vo  may  be  written 
as 

(V^)v0  =  4  Y2  2^v<”c«>cW*  (2-7) 
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To  further  simplify  the  operations,  Equation  2.7  can  be  made  symmetric  with  respect  to  the 
contribution  of  the  vertices  associated  with  each  edge.  This  is  achieved  by  noting  that  if  the 
surface  vectors  of  the  control  volume  sum  to  zero,  the  addition  of  a  constant  (\<t>va)  to  the  above 
formula  does  not  affect  the  gradient  or  divergence.  Thus,  Equation  2.7  may  be  recast  as 

3  1 

(V^)v0  =  -7  -  (4>v,  +  <t>v0)  nco.jCoj*  (2.8) 

Vi, 

Finally,  since  the  area  Am  of  the  median  dual  is  exactly  one  third  that  of  S,  Equation  2.8  becomes: 

(V^)v0  =  1“  ]C  o  +  foo)  nCoiJcojk  (2.9) 

ie{v, 2 

which  is  the  final  formula  used  in  the  reconstruction. 

Minimum-energy  or  "least  squares"  reconstruction  provides  an  alternative  method  for  es¬ 
timating  the  solution  gradient  within  each  cell.  This  reconstruction  process  seeks  to  find  the 
gradient  vector  which  minimizes  the  least  square  error  with  respect  to  the  integral  cell  averages 
of  the  distance  one  neighbors.  Hie  procedure  involves  solution  of  a  (usually)  overconstrained 
system  of  linear  equations.  The  algorithm  can  be  implemented  on  an  edge-by-edge  basis  at  a  cost 
comparable  to  that  of  the  Green-Gauss  formulation.  Details  of  the  least  squares  procedure  may 
be  found  in  Refs.  [9],  [14],  [15]  and  [16],  As  with  the  Green-Gauss  gradient  estimation,  the  set 
of  support  vertices  for  any  node  includes  all  distance  one  neighbors  of  that  node.  However,  the 
least  squares  process  generally  deemphasizes  more  distant  data  as  compared  to  the  Green-Gauss 
method. 

Figure  2  contains  a  quadrilateral,  right  triangular,  and  a  (nearly)  equilateral  triangular  mesh. 
This  figure  displays  not  only  the  physical  mesh,  but  also  the  median  duals  of  each  tessellation. 
These  three  meshes  will  be  a  basis  for  many  of  the  investigations  which  follow.  In  considering 
the  application  of  Equation  2.9  to  unstructured  meshes  containing  polygons  other  than  triangles, 
consider,  for  example,  its  evaluation  on  quadrilateral  cells  depicted  at  the  top  of  Figure  2, 

(V0Ho  ~  ^  2  ^v,nc0,nCoi23  +  ^v3nc0,aCo}«5+ 
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Figure  2:  Physical  mesh  and  median  dual  for  quadrilateral,  right  triangular  and  equilateral 
triangular  tessellations  showing  support  for  gradient  calculation  at  control  volume 
surrounding  node  V0. 


^Vs^CcusCoser  +  ^VV^Coj^Com} 


(2.10) 


where  the  contribution  of  4>v0  has  been  explicitly  cancelled.  Since  the  relationship  in  Equation  2.S 
is  unique  to  triangles,  this  formula  corresponds  only  to  a  midpoint  evaluation  of  the  line  integral 
on  the  dotted  path  shown  in  Figure  2.  As  a  result,  some  of  the  favorable  properties  associated 
with  the  trapezoidal  integration  formula  have  been  sacrificed  and  slightly  less  tolerance  to  mesh 
irregularities  may  be  expected  when  using  the  midpoint  formula  on  quadrilateral  control  volumes. 

Once  the  gradients  of  the  conserved  quantities  are  known,  the  scalars  may  be  reconstructed 
throughout  the  cell: 

<t>(x,y)  =  <j>v0  +  V<f>'r  (2.11) 

where  r  is  a  general  vector  pointing  from  Vo  to  a  point  (x,  y).  However,  such  a  procedure  will 
not  always  yield  monotonic  behavior.  Consequently,  a  limiter  V  is  employed  to  reduce  the  slope 
where  necessary. 

4>(x,y)  =  <f>v,+V  V<f>\j-r  0<V<1  (2.12) 

Taking  V  identically  equal  to  zero  degenerates  the  algorithm  to  first  order. 

2  J2  Limiters 

Barth’s  [1]  original  slope  limiter  is  a  scalar  computed  by  considering  all  edges  incident  upon  V0 
and  applied  directly  to  Let  be  the  reconstructed  value  without  limiting  at  the  median 
point  Vm  of  edge  V0Vj  using  the  gradient  at  Vo. 

<i>r,m  -  <f>V0  +  •  rv0vm  (2.13) 
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Defining  <f>max  =  max(^y0, 4>vs),  ^m.n  =  min(^v0, 4>v}),  At  ,max  ~  ^mat  -<fao>  Aj  ,mtn  ®  4*min  -<Pvc 
A2  =  <t>T,m  —  4>v0,  the  scalar  v|r0j  associated  with  the  gradient  at  Vo  due  to  edge  VoVj  is 

if  A2  >  0 
min(l,^g*)  if  A2  <  0 
1  if  A2  =  0 


(2.14) 


V  is  constructed  from  the  various  Voj  values:  W  =  min(Voj)  where  j  scans  all  vertices  joined 
to  Vo-  This  method  thus  reduces  the  gradient  in  all  directions  equally  and  pushes  Equation  2.12 
toward  first  order  in  all  spatial  directions.  Such  “scalar”  limiting  reduces  both  the  normal  and 
tangential  components  of  the  gradient  and  results  in  a  relatively  dissipative  process. 

This  observation  suggests  that  a  less  dissipative  limiter  may  be  constructed  by  limiting  only 
the  component  of  v<£  along  the  surface  vector  associated  with  each  edge.  A  “directional”  limiter 
may  be  defined  in  which  W  does  not  act  isotropically.  We  begin  by  resolving  the  gradient  into 
components  normal  (n)  and  tangential  (f )  to  an  edge  of  the  centroid  dual. 


V  <l>  =  (V^  •  n)n  +  (V^  •  t)t  (2.15) 

4 1  is  computed  in  a  manner  similar  to  that  described  above  but  is  now  applied  only  to  the  component 
of  the  gradient  normal  to  the  edge  of  the  dual. 


V<^  =  \l;((V^n)n)  +  (v^-T)f  (2.16) 

The  limited  gradient  still  ensures  a  new  maximum  is  not  created  along  the  edge  and  yet  it  now 
avoids  unnecessarily  degrading  the  tangential  component  of  the  slope. 

In  an  attempt  to  further  reduce  the  severity  of  limiting,  a  face-based  implementation  of  this 
directional  limiting  was  examined.  With  this  approach,  the  limiting  is  performed  locally  for  each 
face  separately  without  influencing  the  gradient  stored  at  the  nodes.  Thus,  the  limiting  at  any  one 
face  of  a  control  volume  on  the  dual  mesh  occurs  independently  of  any  limiting  which  may  be 
necessary  on  the  other  faces.  This  procedure  is  analogous  to  that  on  a  structured  mesh  where  the 
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limiter  is  applied  sequentially  in  the  mesh  directions.  The  process  does  not  preserve  the  mean  of 
the  function  value  within  each  cell,  but  can  dramatically  diminish  the  number  of  flux  calculations 
which  invoke  the  limiter.  This  face  based  approach  may  increase  storage  requirements  somewhat 
since  the  limiter  must  be  computed  on  an  edge-by-edge  basis.  Results  with  this  implementation 
are  unsatisfactory  when  used  in  conjunction  with  Barth’s  original  limiting  procedure,  and  alt 
monotone  solutions  are  produced  with  low  absolute  error  in  one  dimension,  multidimensi  ^al 
results  are  erratic.  “Face-based”  limiting  was  not  investigated  in  conjunction  with  other  limiting 
procedures  and  is  not  considered  further  within  this  work. 

Venkatakrishnan  [7]  recently  proposed  a  new  limiter  designed  specifically  to  enhance  the 
convergence  properties  of  the  base  scheme.  This  increase  in  convergence  comes  at  the  expense 
of  strict  monotone  behavior,  as  the  new  (smooth)  limiter  permits  small  local  overshoots  or 
undershoots  in  the  discrete  solution.  This  limiter  may  be  expressed  as: 


1  (^?.min**2)  Az+ZAIAi.min 

^2  <S|~™52|+i5[]min2i2+<" 


if  A2  >  0 
if  A2  <  0 


(2.17) 


Here,  e  is  a  small  number  to  prevent  division  by  zero  while  e  is  chosen  to  be  a  variable  that 
controls  the  degree  of  limiting  and  depends  on  some  estimate  of  mesh  scale  (c2  =  (if  Ax)3).  Ax 
in  this  formula  is  defined  locally  as  the  diameter  of  the  largest  circle  which  may  be  inscribed  into 
a  particular  control  volume. 


23  Flux  Quadrature 

The  divergence  operator  required  for  the  flux  quadrature  of  Equation  2.1  is  closely  related  to  the 
gradient  and  consequently  follows  a  similar  formulation  for  each  component  of  the  equation.  The 
fluxes  are  computed  in  the  present  work  with  Roe’s  flux-difference  split  method  at  the  midpoint 
of  each  edge.  First,  the  state  of  the  flow  on  either  side  of  the  midpoint  of  each  edge  is  constructed 
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from  the  known  values  of  the  state  at  each  node  and  the  limited  gradients  computed  with  the  above 
procedure.  The  flux  function  is  then  evaluated  and  appropriately  scattered  and  accumulated  to 
form  the  solution  change  at  each  node. 


3.  Results  and  Investigations 


The  investigations  examine  issues  of  accuracy,  convergence  and  efficiency  with  respect  to  variation 
in  the  limiter,  type  of  gradient  estimation  and  the  degree  and  quality  of  the  control  volumes  in 
the  discretized  domain.  These  studies  are  conducted  against  a  backdrop  of  widely  computed 
numerical  test  cases  and  closed  form  analytic  solutions  to  the  governing  equations.  In  an  effort  to 
ensure  that  the  conclusions  drawn  are  general  in  nature,  the  numerical  examples  consider  a  variety 
of  smooth  and  nonsmooth  internal  and  external  flows.  Wherever  possible,  efforts  are  made  to 
compare  die  unstructured  results  with  those  from  a  structured  Roe/MUSCL  solver. 


3.1  Convergence  and  Efficiency 

Supersonic  internal  flow  through  a  channel  with  a  4%  circular  arc  bump  provides  an  initial 
assessment  of  the  convergence  properties  of  the  various  methods.  Convergence  behavior  shows 
very  little  dependence  upon  the  reconstruction  procedure  and  results  are  only  presented  using 
Green-Gauss. 

Figure  3  introduces  this  example  through  density  contours  showing  the  three  tessellations  of  a 
65  x  17  mesh.  This  widely  computed  test  case  considers  Mach  1.4  freestream  flow  which  enters 
the  duct  and  sets  up  an  inviscid  shock  reflection  pattern  within  the  domain  [17].  The  meshes 
shown  in  Figure  3  consist  of  quadrilateral,  right  triangular,  and  equilateral  triangular  elements, 
and  the  cases  were  computed  with  both  1st  order  (V^>  =  0)  upwinding  and  linear  reconstruction 
with  all  three  of  the  limiters  outlined  in  Section  2.2. 

All  the  cases  were  run  at  a  CFL  of  0.75  with  local  time-stepping,  and  on  all  the  meshes,  the 
diffusive  character  of  1st  order  upwind  solutions  is  apparent.  While  all  of  the  methods  capture 
die  overall  lambda  shock  structure  within  the  duct,  the  resolution  of  the  interaction  region  near 
die  trailing  edge  permits  discrimination  among  the  higher  order  discrete  solutions.  On  each  of 
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Figure  3:  Density  contours  for  supersonic,  inviscid  flow  through  channel  with  4%  circular 
arc.  Moo  =  1.4,65  x  17  mesh. 
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the  meshes,  for  example,  the  directional  application  of  Barth's  limiter  improves  the  resolution  of 
this  region  considerably  over  the  scalar  implementation.  This  is  consistent  with  the  observation 
alluded  to  in  Section  2.2,  where  it  was  suggested  that  limiting  only  die  component  normal  to  die 
control  volume  face  would  result  in  a  less  dissipative  operator. 

Figure  4  shows  convergence  histories  in  the  L\  norm  for  this  case  on  all  three  meshes.  All 
of  the  examples  with  Venkatakrishnan’s  limiter  used  a  value  of  the  free  parameter  ( K )  of  2.0. 
On  both  triangular  meshes,  variation  of  this  parameter  by  more  than  about  ±0.5  resulted  in  poor 
convergence.  Convergence  behavior  on  the  quadrilateral  meshes  did  not  begin  to  degenerate  until 
this  parameter  was  increased  beyond  5.0. 

3.1.1  Convergence 

Looking  strictly  at  the  convergence  histories  it  would  at  first  appear  that  neither  Barth’s  original 
limiter  nor  the  directional  limiter  allow  the  calculations  to  converge  convincingly.  However, 
an  examination  of  the  contour  plots  at  several  points  in  the  convergence  history  revealed  that 
these  fluctuations  manifest  themselves  only  as  small  amplitude  wiggling  of  the  contours  about 
a  steady  state  as  has  been  reported  by  other  investigations  [7].  For  the  directional  limiter,  the 
fluctuations  are  about  an  order  of  magnitude  smaller  and  perturbations  of  the  contour  lines  are  not 
immediately  evident  (although  undoubtedly  there).  This  is  consistent  with  the  results  displayed  by 
the  convergence  histories.  The  directional  limiter  enhanced  convergence  by  roughly  an  additional 
order  of  magnitude,  obviously  indicating  less  activity. 

One  plausible  explanation  for  this  behavior  views  these  fluctuations  in  the  convergence  his¬ 
tories  as  the  results  of  sporadic  firing  of  the  limiter,  thus  preventing  absolute  convergence. 
This  assertion  is  further  supported  by  the  work  of  Ref.  [7],  where  the  free  parameter  (K)  in 
Venkatakrishnan’s  limiter  was  introduced  to  prevent  the  limiter  from  prematurely  reacting  to 
slight  oscillations.  In  the  present  examples,  incorporation  of  this  tolerance  facilitates  convergence 
to  machine  zero  on  all  three  meshes  (32-bit  machine).  The  directional  limiter  displays  behavior 


Figure  4:  Convergence  histories  for  various  schemes  on  each  mesh  for  supersonic  inviscid 
flow  through  channel  with  4%  circular  arc.  CFL  =  0.75 
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similar  to  that  of  Barth.  However,  in  all  of  the  cases ,  the  fact  that  it  only  retards  one  component 
of  the  gradient  results  in  somewhat  better  convergence  since  the  limiting  is  less  severe.  This  issue 
will  be  revisited  in  die  next  investigation,  where  the  behavior  of  each  of  die  limiters  is  mapped 
out  in  detail. 

3.U  A  Note  on  Efficiency 

As  a  final  note,  it  is  worth  mentioning  that  while  all  three  meshes  require  roughly  die  same  number 
of  iterations  to  converge,  the  calculations  on  the  triangular  meshes  were  50%  more  expensive  in 
both  memory  and  time  than  those  on  the  quadrilateral  meshes.  This  follows  since  the  scheme 
proceeds  on  an  edge-by-edge  basis  and  six  edges  (instead  of  four)  are  incident  upon  each  vertex 
in  the  interior  of  the  mesh.  In  three  dimensions,  a  typical  interior  vertex  has  six  incident  edges  on 
hexahedral  meshes,  while  a  corresponding  vertex  will  have  approximately  12-14  incident  edges 
on  a  tetrahedral  mesh.  Thus,  about  2  to  2.5  times  times  the  storage  will  be  needed  if  a  tetrahedral 
mesh  is  chosen.  The  hope  is,  of  course,  that  the  additional  edges  and  flux  evaluations  in  triangular 
and  tetrahedral  domains  may  enhance  the  wave  propagation  within  the  discrete  domain  and  lead  to 
more  accurate  numerical  solutions.  Nevertheless,  an  examination  of  Figure  3  does  not  reveal  any 
obvious  benefit  stemming  from  the  additional  edges  present  in  the  tessellation,  and  it  is  difficult 
to  justify  the  additional  expense  based  strictly  upon  accuracy  arguments. 


3.2  Order  of  Accuracy  and  Absolute  Error 

3.2.1  Supersonic  Vortex 

While  the  example  in  the  previous  section  gives  some  basis  for  comparing  the  methods,  quanti¬ 
tative  measurement  of  the  order  of  accuracy  and  discretization  error  associated  with  each  scheme 
is  best  performed  on  a  test  case  for  which  an  exact,  closed  form,  analytic  solution  exists.  Exami- 
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nation  of  a  two  dimensional  supersonic  vortex  provides  just  such  an  opportunity.  Since  this  is  a 
shock  free  compressible  flow,  the  measured  order  of  accuracy  is  not  corrupted  by  limiter  action 
near  shocks  and  die  behavior  will  correspond  to  what  one  may  expect  of  each  method  in  smooth 
regions  of  the  flow.  Using  this  case  as  a  diagnostic  tool  permits  direct  examination  of  the  effects 
of  polygons  of  various  shapes  and  degrees,  and  also  provides  an  opportunity  to  study  the  tolerance 
of  the  methods  to  poor  quality  meshes. 

The  inviscid,  isentropic,  supersonic  flow  of  a  compressible  fluid  between  concentric  circular 
arcs  presents  a  flow  in  which  the  velocity  varies  inversely  with  the  radius.  This  flow  is  a  particularly 
useful  test  case  for  upwind  methods  since  the  numerical  solutions  must  propagate  infinitesimally 
weak  waves  accurately  to  perform  the  large  turning  without  disrupting  the  radial  flow  distribution 
or  introducing  shocks.  This  flow  is  also  of  some  practical  interest  as  it  has  been  used  as  a  segment 
of  the  flow  distribution  in  designing  passages  for  the  supersonic  blading  of  compressors  and 
turbines  (18, 19],  as  well  as  for  a  supersonic  through-flow  fan  stage  [20]. 

The  expression  for  density  p  as  a  function  of  radius  r  is  given  by: 


p(r)  —  Pi 


(3.1) 


where  ilf,  and  r,  are  the  Mach  number  and  the  radius  at  the  inner  arc.  To  evaluate  the  order  of 
accuracy  of  each  of  the  methods,  a  series  of  discrete  solutions  to  this  test  case  is  obtained  on  the 
quadrilateral  as  well  as  on  equilateral  and  right  triangular  meshes.  For  each  tessellation,  solutions 
are  sought  on  a  set  of  three  telescoping  grids  with  31  x  6, 61  x  1 1,  and  121  x  21  nodes.  The  Mach 
number  at  the  inner  radius  r,  is  specified  at  2.25  and  the  outer  radius  r0  at  1 .384  r, .  Figure  5  shows 
the  three  sets  of  regular  meshes  used  in  the  simulation.  For  these  regular  meshes,  the  aspect  ratio 
is  of  order  1. 

The  simulations  are  initiated  by  releasing  the  inlet  profile  into  a  nearly  evacuated  duct  and 
converged  to  steady  state.  Figure  6  displays  pressure  contours  (inc.=  0.25)  resulting  from  the 
computations  on  the  medium  meshes.  The  extreme  diffusivity  of  the  first  order  schemes  again 
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becomes  apparent  in  these  plots,  and  the  pressure  in  the  first  order  simulations  drops  dramatically 
along  stream  lines  as  the  flow  proceeds  around  the  duct.  Some  evidence  of  dissipation  is  also 
evident  in  the  higher  order  discrete  solutions  but  it  is  not  nearly  as  pronounced  as  in  the  first  order 
cases.  Convergence  behavior  for  all  the  schemes  mimics  that  for  the  duct  flow  case  presented 
earlier.  The  noise  in  the  solutions  with  the  first  two  limiters  will  be  discussed  shortly. 

The  existence  of  an  exact  solution  to  this  problem  makes  it  possible  to  directly  compute 
distributions  of  error  throughout  the  domain. 

Figure  7  contains  contour  plots  of  density  error  for  the  medium  triangular  and  quadrilateral 
meshes  for  all  the  methods.  Each  contour  level  indicates  2%  error  in  the  local  density  distribution. 
Note  that  the  error  contours  do  not  display  any  anomalies  approaching  the  walls  or  boundaries, 
suggesting  that  the  accuracy  near  the  wall  is  the  same  as  that  in  the  interior  of  the  field. 

By  comparing  the  error  in  the  discrete  solutions  on  a  successively  refined  sequence  of  tele¬ 
scoping  meshes,  quantitative  measurements  of  both  order  of  accuracy  and  absolute  error  are 
possible.  Note  that  since  the  simulations  were  run  to  convergence,  the  error  measured  is  r  t 
directly  truncation  error  but  rather  the  error  in  the  discrete  solutions.  The  Appendix  tabulates 
the  L\  and  L2  norm  of  the  density  error  in  the  discrete  solutions  on  all  the  meshes  considered  in 
this  example.  In  addition,  it  offers  order  of  accuracy  estimates  for  the  procedures  on  each  type 
of  polygon.  In  general,  the  L\  and  L2  norms  behave  quite  similarly,  and  the  order  of  accuracy 
of  the  discrete  solutions  is  roughly  the  same  in  using  either  measure.  Since  the  L2  norm  is  more 
sensitive  to  extrema,  this  fact  suggests  that  discretization  error  is  being  driven  out  of  the  domain 
uniformly,  and  supports  the  contention  that  the  boundaries  are  free  from  anomalies.  Selected  data 
will  be  extracted  from  these  tables  to  aid  in  the  analysis  of  the  following  sections. 

3.2.2  Polygon  Degree 

The  first  investigation  examines  the  influence  of  the  shape  and  degree  of  the  polygons  in  the 
domain  on  the  accuracy  of  the  discrete  solutions.  Since  results  free  from  the  effects  of  limiters 


21 


are  sought,  the  simulations  are  started  from  the  exact  solution,  and  discrete  solutions  are  obtained 
without  the  use  of  flux  limiters.  The  order  of  accuracy  of  the  discrete  method  is  assessed  by 
performing  a  linear  regression  on  log-log  plots  of  absolute  density  error  in  both  L\  and  L2  norms 
versus  normalized  grid  spacing  represented  by  the  reciprocal  square  root  of  the  number  of  nodes 
in  the  domains.  Figure  8  presents  the  order  of  accuracy  of  the  methods  in  L\,  and  Table  A.1  in 
the  Appendix  contains  the  data  used  to  generate  this  chart. 

In  this  example,  both  the  order  of  accuracy,  and  absolute  error  on  the  quadrilateral  and 
equilateral  triangular  tessellations  are  comparable.  The  order  of  accuracy  on  these  two  meshes  is 
S — 10%  better  than  that  on  the  right  triangular  mesh  with  either  reconstruction  procedure.  The  first 
order  results  show  a  more  pronounced  deficit,  placing  the  order  of  accuracy  of  the  right  triangular 
mesh  at  about  0.7  as  compared  to  1 . 1  and  1 .0  for  die  quads  and  equilateral  triangles,  respectively. 
This  fact  suggests  that  something  in  the  very  nature  of  the  right  triangles  is  responsible  for  the 
lower  order  of  accuracy  on  such  meshes. 

The  sketches  in  Figure  9  permit  a  direct  comparison  of  the  dual  meshes  on  quadrilateral  and 
right  triangular  tessellations.  In  examining  this  figure  it  becomes  apparent  that  the  face  a,  which 
is  pierced  by  the  edge  LR  is  at  a  right  angle  a  on  the  quadrilateral  mesh  but  not  on  the  stretched 
triangular  mesh.  (Note  that  Figure  2  shows  that  a  is  a  right  angle  on  equilateral  triangular  meshes 
as  well.)  In  fact,  as  the  mesh  stretching  increases,  the  alignment  of  faces  like  a,  b ,  d,  and  e  with 
their  respective  edges  becomes  less  and  less  orthogonal. 

When  the  first  order  essentially  1-D  scheme  evaluates  the  flux  across  a  poorly  aligned  face, 
like  a,  the  Riemann  solver  expects  data  normal  to  the  face,  but  the  data  provided  by  either  vertex 
of  edge  LR  is  far  removed  from  the  true  normal,  and  necessarily  introduces  an  error  into  the  flux 
evaluation.  On  unskewed  quadrilateral  meshes,  or  equilateral  triangular  meshes  a  is  very  nearly  a 
right  angle.  As  a  result,  the  data  introduced  into  the  approximate  Riemann  solver  is  well  aligned 
with  the  normal  to  the  face,  and  such  an  error  is  avoided.  Such  arguments  support  the  hypothesis 
that  the  first  order  scheme  may  degrade  rapidly  on  highly  stretched  right  triangular  meshes. 
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Order  of  Accuracy  in  LI 
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Regular  Rt.  Equilat. 

Quads  Triangles  Triangles 


ID  1st  Order 
I  GG  -  No  Limiter 
E3  L2  -  No  Limiter 


Figure  8:  Order  of  accuracy  of  unlimited  schemes  on  regular  meshes  for  the  supersonic  vortex 
problem 
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States  Used  After  Reconstruction 

Figure  9:  Setting  up  of  the  local  Riemann  problem  on  a  face  of  the  auxiliary  meshes  for 
quadrilateral  and  right  triangular  tessellations 
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When  linear  reconstruction  is  introduced  into  the  formation  of  the  numerical  flux,  die  new  left 
and  right  states  for  the  Riemann  problem  are  now  formed  on  either  side  of  face  a  in  Figure  9.  These 
states  are  labeled  L'  and  R  in  the  sketch  and  take  into  account  linear  variation  within  the  auxiliary 
cell.  An  alternate  view  of  die  situation  is  to  realize  that  the  reconstruction  step  has  taken  data 
into  account  from  all  of  the  distance  one  vertices  adjacent  to  L  and  R,  thus  enlarging  die  stencil 
to  permit  information  to  propagate  in  the  face-normal  direction.  Since  there  is  no  direct  edge 
connection  normal  to  face  a,  the  first  order  scheme  necessarily  makes  an  error  when  propagating 
information  against  the  mesh  diagonal.  Reconstruction  allows  the  higher  order  scheme  to  recover 
more  than  the  single  order  of  magnitude  associated  with  the  slope  estimation  and  very  nearly 
matches  the  discrete  solutions  of  the  quadrilaterals  and  equilateral  triangles.  By  enlarging  the 
stencil  to  incorporate  data  in  the  face-normal  direction,  reconstruction  reduces  the  misalignment 
of  the  Riemann  problem  present  in  the  first  order  discrete  solutions. 

Some  general  comments  stem  from  this  investigation  with  the  unlimited  schemes  on  regular 
meshes.  (1)  On  regular  meshes  the  Green-Gauss  and  least  squares  gradient  estimation  procedures 
yield  similar  results.  (2)  The  regular  quadrilateral  and  equilateral  meshes  yield  very  similar 
discrete  solutions,  both  in  order  of  accuracy,  and  absolute  magnitude  of  the  error.  Nevertheless, 
the  simulations  on  quadrilaterals  require  approximately  50%  less  storage  and  CPU  simply  by 
virtue  of  the  fact  that  fewer  edges  exist  in  the  domain.  (3)  The  scheme  seems  to  treat  right 
triangles  as  distorted  equilateral  triangles,  and  increases  the  reliance  on  the  gradient  estimation 
for  the  production  of  accurate  discrete  solutions. 

3.2.3  Limiter  Behavior 

The  next  set  of  numerical  examples  introduces  the  slope  limiters  into  the  simulations  on  the 
three  sets  of  regular  meshes  previously  presented.  This  investigation  is  designed  to  quantitatively 
compare  the  effects  of  limiting  on  scheme  accuracy  with  each  of  the  three  polygons  types. 

As  mentioned  earlier,  the  limiter  originally  proposed  for  the  scheme  [1]  is  a  scalar  correspond- 


26 


ing  to  the  severest  requirement  which  maintains  monotonicity.  The  directional  implementation 
of  this  limiter  was  presented  in  Section  2.2  ami  is  designed  to  reduce  the  severity  of  die  limiting 
by  decreasing  only  die  normal  component  of  the  gradient  The  third  limiter  examined  was  re¬ 
cently  proposed  by  Venkatakrishnan  [7].  This  smooth  limiter  is  less  severe,  and  is  invoked  less 
frequently  than  Barth’s  original  limiter. 

After  converging  numerical  simulations  on  the  medium  quadrilateral  mesh  for  about  10, 000 
iterations  to  reach  a  steady  state,  the  behavior  of  the  limiter  with  each  scheme  was  monitored 
for  several  hundred  iterations.  Table  1  contains  an  excerpt  from  this  data  and  tracks  the  mean 
behavior  of  each  limiter  for  several  iterations.  The  first  two  columns  display  the  percentage  of 
all  edges  within  the  domain  on  which  the  limiter  was  applied  for  each  reconstruction  method. 
The  last  two  columns  tabulate  the  average  value  of  the  limiter  on  only  those  faces  where  the 
limiter  was  active,  thus  giving  a  picture  of  the  severity  of  the  limiting  taking  place.  Of  course, 
since  the  discrete  answers  are  slightly  different,  the  edges  with  flux  evaluations  that  invoke  the 
limiter  are  not  always  the  same  from  solution  to  solution.  The  lack  of  absolute  convergence 
with  Barth’s  limiter  is  immediately  apparent  in  the  first  two  columns  as  the  number  of  limited 
edges  changes  from  time  step  to  time  step  but  hovers  around  17%.  The  average  value  of  the 
limiter  on  these  edges  is  around  0.25  indicating  a  very  nearly  first  order  flux  evaluation  on  these 
edges.  Applying  this  limiter  only  to  the  normal  component  of  the  gradient  substantially  reduces 
the  magnitude  of  the  slope  degradation  as  shown  by  the  tabulated  results  with  the  directional 
limiter.  With  this  vector  application  of  Barth’s  limiter,  the  limited  slopes  retain  92  -  94%  of  their 
magnitude  and  monotonicity  is  still  guaranteed.  The  fact  that  the  exact  number  of  limited  edges 
is  now  substantially  more  stable  gives  evidence  of  the  deeper  convergence  afforded  by  the  milder 
limiting.  Venkatakrishnan’s  limiter  is  designed  to  fire  less  frequently,  and  this  table  shows  that 
with  K  =  10  it  is  triggered  on  only  5  -  6%  of  the  edges  depending  on  the  reconstruction,  and 
reduces  die  slopes  by  an  average  of  only  2  -  3%. 

Figure  10  summarizes  the  effects  of  limiting  on  the  order  of  accuracy  with  the  regular 
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Table  1:  Table  showing  severity  and  frequency  of  limiter  filing  after  practical  convergence 
average  values  on  medium  mesh  for  supersonic  vortex  flow 

Barth’*  limiter 


itt  . . . :■ 

Mean  W 

■  —  ■ 

GG  LS  j 

1832677  16.62402 

17.01772  15.99409 
16.63386  18.34646 
19.19291  16.70275 

17.70669  15.86614 
16.80118  17.42126 
1838268  15.75787 

17.14567  1539055 
16.68307  17.02756 
1834803  15.93504 

03444862  03487914 
03593847  0.2473031 
03549107  0.2249657 
03311892  03505058 
03648712  03523419 
0.2390375  03494676 
0.2293676  03628614 
03513207  03808488 
0.2457611  03357707 
0.2329683  03580887 

Directional  limiter 


%  edges  limited 

Mean  V 

GG  LS 

GG  LS 

24.72441  24.58661 

24.72441  24.57677 

24.72441  24.57677 

24.72441  24.57677 

24.72441  24.55709 

24.72441  24.56693 
24.72441  24.56693 

24.72441  24.56693 

24.72441  24.54725 
24.72441  24.56693 

0.9240345  0.9457827 
0.9240335  0.9457660 
0.9240359  0.9457663 
0.9240357  0.9457676 
0.9240360  0.9457363 
0.9240357  0.9457521 
0.9240364  0.9457527 
0.9240357  0.9457538 
0.9240348  0.9457356 
0.9240317  0.9457530 

Venkatakrishnan  limiter,  K  =  10 


%  edges  limited 

Mean  V 

GG  LS 

GG  LS  | 

6.072834  4.970472 
6.072834  4.970472 
6.072834  4.970472 
6.072834  4.970472 
6.072834  4.970472 
6.072834  4.970472 
6.072834  4.970472 
6.072834  4.970472 
6.072834  4.970472 
6.072834  4.970472 

0.9779245  0.9867367 
0.9779246  0.9867366 
0.9779246  0.9867365 
0.9779242  0.9867367 
0.9779242  0.9867365 
0.9779242  0.9867368 
0.9779245  0.9867368 
0.9779246  0.9867369 
0.9779245  0.9867367 
0.9779245  0.9867367 
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Figure  10:  Order  of  accuracy  of  limited  schemes  on  regular  polygons  for  the  supersonic  vortex 
problem 
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meshes  of  each  polygon  shape.  These  results  summarize  the  data  contained  in  Table  A.2  in  die 
Appendix  which  presents  both  error  norms  and  order  of  accuracy  estimates.  Figure  10  shows 
that  Venkatakrishnan’s  limiter  nearly  reproduces  the  behavior  of  the  unlimited  schemes  for  this 
smooth  flow.  This  is  consistent  with  the  high  values  of  the  limiter  observed  in  Table  1. 

In  cases  with  no  limiting,  or  with  very  little  limiting  -  as  with  Venkatakrishnan’s  limiter  -  the 
right  triangular  elements  almost  match  the  accuracy  of  the  other  two  tessellations.  However,  this 
has  come  from  an  increased  reliance  on  the  gradient  estimation.  The  behavior  with  die  directional 
limiter  provides  a  good  demonstration  of  this.  When  the  gradient  at  L  or  R  is  decreased  due  to 
limiting,  the  stencil  will  again  revert  to  the  poorly  aligned  first  order  Riemann  solution  and  the 
propagation  of  information  will  be  restricted.  Thus,  the  reconstructed  values  at  the  cell  face  will 
degrade  rapidly.  In  Equation  2.13,  accurate  estimations  of  <f>r,m  rely  ever  more  heavily  on  \j<f> 
since  4>v0  represents  data  not  normal  to  the  face.  Thus,  limiting  xj<t>  will  degrade  the  discrete 
solution  more  severely  than  on  well  aligned  meshes.  On  the  finest  right  triangular  mesh  with  the 
directional  limiter  the  error  is  nearly  10  times  that  on  the  equilateral  triangles  or  quadrilaterals 
with  the  same  limiter  (see  Table  A.2  in  Appendix). 

32.4  Polygon  Quality 

While  it  is  informative  to  evaluate  the  scheme  performance  on  the  meshes  in  the  previous  inves¬ 
tigation,  such  smooth,  regular,  aspect  ratio  ~  1  polygons  are  rarely  found  in  practice.  Having 
established  a  reference  level  for  scheme  accuracy,  focus  now  shifts  to  discrete  solutions  on 
stretched  and  distorted  meshes. 

Figure  11  displays  two  sets  of  telescoping  stretched  meshes.  The  quadrilateral  and  triangular 
elements  were  formed  on  the  same  set  of  vertices  and  the  aspect  ratio  of  the  (quadrilateral) 
elements  is  about  40.  The  original  scalar  limiter,  and  the  directional  limiter  refused  to  converge 
on  the  coarse  and  medium  triangular  meshes  for  this  case,  so  Table  A.2  in  the  Appendix  only 
contains  results  with  the  first  order  scheme  and  Venkatakrishnan’s  limiter. 
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Order  of  Accuracy  in  LI 


1st  Order  GG-  LS-  Structured 
Venkata*  Venkata*  Roe 
krishnan  krishnan 


Figure  12:  Order  of  accuracy  of  limited  schemes  on  stretched  polygons  for  the  supersonic 
vortex  problem 
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Figure  12  contains  a  comparison  of  the  order  of  accuracy  information  derived  from  die 
tabulated  data.  The  first  order  results  place  the  order  of  accuracy  on  the  quadrilateral  cells  at 
0.93  and  the  right  triangular  scheme  at  only  0.41.  Such  results  are  consistent  with  the  argument 
presented  earlier  which  predicted  that  the  alignment  problem  becomes  worse  with  increasing  cell 
aspect  ratio. 

At  first  glance,  the  improvement  in  order  of  accuracy  for  the  reconstructed  solutions  on  the 
triangular  grid  is  very  impressive.  Although  the  order  of  accuracy  is  still  5  -  10%  lower  than 
on  the  quadrilateral  meshes,  reconstruction  has  improved  the  solution  by  more  than  a  full  order 
of  magnitude.  However,  after  examining  the  data  in  Table  A.3,  it  is  clear  that  even  on  the  finest 
mesh,  the  discrete  solutions  on  the  right  triangular  mesh  exhibit  4-10  times  more  error  than  on 
the  quadrilateral  tessellation  of  the  same  set  of  nodes.  Despite  this,  the  discrete  solutions  on  the 
fine  triangular  mesh  are  quite  reasonable,  and  on  a  sufficiently  fine  mesh,  the  gradient  estimation 
is  able  to  compensate  for  the  error  in  the  first  order  scheme. 

The  right  triangle’s  increased  dependency  on  the  accuracy  of  the  gradient  estimation  becomes 
apparent  when  studying  the  behavior  of  the  discrete  solutions  on  the  coarse  meshes.  On  these 
meshes,  only  5  points  spanned  the  domain  in  the  stream  wise  direction.  The  resulting  gradient 
estimates  led  to  discrete  solutions  which  are  actually  worse  than  the  first  order  simulations.  On 
the  quadrilateral  meshes,  similar  gradient  estimates  improved  the  discrete  solutions  by  a  factor  of 
2-6. 

Finally,  in  regarding  Table  A.3,  note  that  the  discrete  solution  on  the  right  triangular  mesh  using 
least  squares  reconstruction  converged  to  a  result  which  unstarted  the  flow.  Since  it  converged 
to  a  different  physical  solution,  this  point  was  not  considered  in  the  slope  estimate  provided 
in  the  table.  The  results  for  this  test  case  appear  to  substantiate  observations  made  previously. 
However,  the  schemes  poor  performance  on  the  coarse  right  triangular  mesh  suggests  that  they 
should  be  regarded  as  preliminary  until  the  results  can  be  substantiated  by  computations  on  still 
finer  stretched  meshes. 
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Estimated  Order  of  A  -  yin  LI 


Figure  14:  Order  of  accuracy  of  limited  schemes  on  randomly  distorted  polygons  for  the  su¬ 
personic  vortex  problem 
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In  order  to  examine  the  schemes  tolerance  to  distorted  elements,  the  regular  meshes  shown 
in  Figure  5  were  perturbed  to  introduce  localized  skewing.  Figure  13  displays  the  resulting 
distorted  meshes.  Each  mesh  point  was  randomly  displaced  within  a  small  region  of  its  local 
neighborhood.  The  resultant  sets  of  nodes  were  then  connected  to  form  both  quadrilateral  and 
triangular  polygons.  Notice  that  these  meshes  do  not  telescope,  since  the  magnitude  of  the  mesh 
point  displacement  scales  with  the  normalized  mesh  spacing. 

Figure  14  compares  the  estimated  order  of  accuracy  for  each  discrete  method  on  the  two  poly¬ 
gon  types.  This  chart  is  again  drawn  from  the  information  provided  in  the  Appendix  (Table  A.4). 
The  chart  displays  a  clear  degradation  in  accuracy  on  the  quadrilateral  meshes.  The  triangular 
elements,  in  contrast,  produce  results  which  nearly  match  those  on  regular  equilateral  triangles 
presented  earlier.  Table  A.4  shows  that,  in  general,  somewhat  higher  discretization  error  exists 
on  the  distorted  triangular  meshes.  However,  this  table  also  shows  that  the  least  squares  cases  are 
extremely  tolerant  to  the  mesh  distortion  and  the  error  in  these  solutions  is  very  nearly  as  low  as 
that  on  the  regular  triangular  meshes. 

These  observations  support  several  general  statements  about  the  method’s  performance  on 
distorted  polygons.  It  appears  that  on  both  tessellations,  the  least  squares  gradient  estimation  is 
far  more  capable  of  producing  reliable  results  on  distorted  meshes.  Moreover,  its  also  evident 
that  the  trapezoidal  integration  of  the  Galerkin  portion  of  the  numerical  flux  function  on  triangles 
makes  these  polygons  far  more  tolerant  of  mesh  distortion.  The  order  of  accuracy  on  the  distorted 
quadrilateral  meshes  degenerates  much  more  quickly,  since  the  edge  formulas  only  result  in 
midpoint  quadrature  for  these  polygons.  Notice  however,  that  the  absolute  magnitude  of  the  error 
on  all  of  the  quad  meshes  remains  very  low  which  suggests  that  the  coefficient  of  the  discretization 
error  expression  remains  small. 
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33  External  Inviscid  Flow  -  Transonic  NAC  A  0012  Airfoil 

Transonic,  inviscid  flow  past  a  NACA  0012  airfoil  poses  a  challenging  and  realistic  problem  on 
which  to  verify  the  features  of  the  various  algorithms  discussed  in  the  preceding  sections.  Mach 
0.8  flow  at  an  angle  of  attack  of  1.25*  results  in  a  steady  solution  with  a  slip  line  at  the  trailing 
edge,  a  strong  shock  on  the  suction  surface,  and  a  weak  shock  on  the  pressure  side.  The  upper  and 
lower  frames  on  the  left  of  Figure  IS  show  the  quadrilateral  and  triangular  meshes  (respectively) 
used  in  the  simulation.  The  remainder  of  this  figure  contains  a  sampling  of  typical  results  for 
several  combinations  of  limiter  and  reconstruction  method.  Results  on  the  quadrilateral  control 
volumes  are  displayed  above  the  discrete  solutions  on  the  triangular  meshes. 

The  mesh  employed  in  this  example  consists  of  131  x  65  vertices,  and  the  circumferential 
resolution  was  chosen  such  that  the  lower  shock  just  formed  when  using  the  Green-Gauss  recon¬ 
struction  in  conjunction  with  Venkatakrishnan’s  limiter,  on  the  quadrilateral  mesh.  The  second 
column  in  Figure  15  contains  the  first  order  discrete  solutions  on  each  mesh.  As  in  the  previous 
examples,  the  first  order  triangular  scheme  produces  a  result  which  appears  more  diffuse  than  that 
on  the  quad  mesh.  The  triangles  near  the  surface  are  very  nearly  right  triangles  with  an  aspect 
ratio  of  roughly  10,  and  the  most  likely  explanation  is  again  poor  alignment  of  the  Riemann  solver. 

The  last  2  columns  of  Figure  15  display  reconstructed  solutions  using  the  directional  limiter 
with  Green-Gauss  and  Venkatakrishnan’s  limiter  with  least  squares.  Reconstruction  dramatically 
improves  the  quality  of  the  discrete  solutions  not  only  because  of  the  higher  order  of  accuracy 
of  the  resulting  scheme,  but  also  as  a  result  of  the  extended  support  stencil  used  in  forming  the 
Riemann  problem.  Ultimately,  the  discrete  solutions  on  the  two  meshes  are  extremely  close. 
In  fact,  with  Venkatakrishnan’s  limiter  and  least  squares  reconstruction,  the  solutions  with  the 
quadrilaterals  and  the  triangles  are  very  nearly  indistinguishable  and  the  lift  coefficients  match  to 
3  significant  digits. 

The  Cp  plots  in  Figure  16  provide  a  more  quantitative  assessment  of  the  different  schemes. 
This  figure  contains  the  results  from  all  tested  combinations  of  limiter  and  reconstruction  method. 
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Figure  15:  Meshes  and  isobars  for  a  sampling  of  typical  discrete  solutions  for  a  NACA  0012  at 
-  0.8  and  a  =  1.25  deg.  (upper -quadrilaterals,  lower -triangles,  131  x  65 
mesh) 


Also  included  for  reference  are  results  from  a  cell-centered  structured  Roe/MUSCL  solver  [21]. 
The  profiles  are  strongly  dependent  upon  the  capturing  of  the  weak  lower  shock.  Less  accurate 
methods  fail  to  predict  this  shock  and  represent  it  as  simply  a  smooth  compression.  With  die 
exception  of  the  structured  mesh  results,  hollow  squares  denote  data  with  the  quadrilateral  meshes, 
while  filled  triangles  represent  the  triangular  elements. 

In  the  first  order  solutions,  only  the  quadrilateral  mesh  gives  any  indication  of  the  lower 
shock’s  existence.  On  either  of  these  meshes,  the  original  scalar  limiter  fails  to  predict  the  shock 
as  well.  Indeed,  these  results  did  not  converge  convincingly,  and  the  residual  stalled  after  dropping 
approximately  3  orders  of  magnitude.  The  directional  implementation  of  this  limiter  reduced  the 
residual  by  an  additional  two  orders,  and  the  lower  shock  is  evident  in  the  discrete  solutions  on 
both  tessellations.  As  in  the  previous  examples,  Venkatakrishnan’s  limiter  converges  to  machine 
zero,  but  over/undershoots  appear  in  all  the  examples. 

The  behavior  of  the  two  reconstruction  methods  follows  the  same  trends  as  in  the  shockless 
flows  considered  earlier.  Least  squares  appears  to  produce  more  accurate  gradient  estimates,  and 
the  advantage  is  particularly  evident  on  the  triangular  meshes.  This  method  also  shows  a  reduced 
sensitivity  to  polygon  shape  because  it  more  appropriately  weights  the  data  surrounding  a  vertex. 
On  triangular  meshes  formed  by  division  of  quadrilaterals,  the  Green-Gauss  reconstruction  may 
actually  introduce  a  diagonal  bias,  while  the  least  squares  reduces  the  weighting  of  these  points 
due  to  their  distance  from  the  central  vertex.  On  none  of  these  relatively  smooth  grids  did 
the  additional  edge  calculations  present  on  the  right  triangular  meshes  resul  n  any  apparent 
advantage. 
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Figure  16:  CP  distributions  with  various  schemes  for  Mach  0.8,  a  =  1.25°  NACA  0012  airfoil 
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4.  Conclusions 


The  accuracy  and  efficiency  of  a  variety  of  common  edge-based  reconstruction  schemes  are 
examined  on  unstructured  meshes  with  various  limiters  and  types  of  polygons.  For  meshes  with 
cell  aspect  ratios  near  one,  the  accuracy  of  the  discrete  solutions  on  triangles  and  quadrilaterals 
is  nearly  the  same  in  both  order  and  absolute  magnitude.  Right  triangular  elements,  however, 
appear  to  be  viewed  as  distorted  equilateral  triangles  -  with  an  associated  reduction  in  scheme 
accuracy  which  is  particularly  severe  in  the  piecewise  constant,  first  order  scheme.  The  situation 
is  exacerbated  when  the  mesh  is  stretched  to  produce  high  aspect  ratio  elements  and  the  Riemann 
solver  is  more  poorly  aligned  with  the  edge  which  introduces  the  first  order  data.  Higher-order 
calculations  on  such  meshes  rely  on  the  reconstruction  to  extend  the  support  of  the  first  order 
stencil.  This  results  in  an  increased  burden  on  the  gradient  estimation,  and  makes  them  easily 
degraded  by  limiting.  Even  when  the  reconstruction  helps  recover  the  loss  in  order  of  accuracy, 
the  absolute  level  of  error  in  the  discrete  solutions  remains  a  factor  of  5-10  higher  than  on  the 
quadrilateral  meshes.  The  fact  that  quadrilaterals  are  not  rigid  figures  makes  it  possible  to  stretch 
these  cells  without  introducing  skewing,  and  elevated  cell  aspect  ratio  does  not  degrade  these 
discrete  solutions.  Nevertheless,  on  poor  quality  meshes,  the  midpoint  integration  and  smaller 
support  stencil  of  the  quadrilateral  meshes  combine  to  make  the  discrete  solutions  degenerate 
more  rapidly. 

A  new  directional  implementation  of  Barth’s  original  scalar  limiter  is  introduced  which  pre¬ 
serves  the  cell  average  and  retains  monotonicity.  This  approach  limits  only  the  component  of 
the  gradient  normal  to  the  face.  Numerical  experiments  show  that  this  implementation  signif¬ 
icantly  reduces  the  dissipation  introduced  by  the  original  scalar  procedure.  Venkatakrishnan’s 
smooth  limiter  performs  as  expected  and  converges  for  all  test  cases  while  producing  small 
over/undershoots  in  the  discrete  solutions.  A  very  promising  approach  toward  further  reducing 
the  introduction  of  excessive  dissipation  into  the  solutions  may  be  the  directional  implementation 
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of  Venkatakri  shnan  ’  s  limiter. 

On  regular  and  stretched  meshes,  the  additional  edges  in  the  triangular  meshes  do  not  lead  to 
any  apparent  accuracy  advantage  over  the  quadrilateral  scheme.  These  additional  edges,  however, 
do  mandate  50%  higher  (in  2-D)  storage  and  CPU  requirements  (over  200%  in  3-D)  -  simply  by 
virtue  of  the  tessellation  and  the  fact  that  the  scheme  proceeds  on  an  edge  basis.  These  observations 
support  a  possible  mesh  generation  strategy  which  removes  unnecessary  edges  from  boundary 
layer  regions  and  other  regular  portions  of  an  unstructured  triangular  mesh,  and  processes  such 
regions  as  a  collection  of  mixed  quadrilateral  and  triangular  elements. 
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